pip install MDAnalysis

import MDAnalysis as mda
from MDAnalysis.analysis import rms

# Load the structure and trajectory
u = mda.Universe(r"C:\Users\jakea\Desktop\Molecular Dynamics\v2_pro8_progress\wt_pro8_trajout\traj.pdb", r"C:\Users\jakea\Desktop\Molecular Dynamics\v2_pro8_progress\wt_pro8_trajout\traj.dcd")

# Select the protein
protein = u.select_atoms("protein")

# Set a reference structure (e.g., first frame)
reference = u.select_atoms("protein").positions

# Prepare a dictionary to store RMSD for each residue
residue_rmsd = {residue.resid: [] for residue in protein.residues}

# Loop over the trajectory
for ts in u.trajectory:
    # Calculate RMSD for each residue
    for residue in protein.residues:
        residue_selection = residue.atoms
        rmsd_value = rms.rmsd(residue_selection.positions, reference[residue_selection.indices])
        residue_rmsd[residue.resid].append(rmsd_value)

# Now residue_rmsd contains RMSD values for each residue across the trajectory

----------------------------------------------------------------------------------------------------------------


import csv

# Write RMSD values to a CSV file
with open(r"C:\Users\jakea\Desktop\Molecular Dynamics\v2_pro8_progress\wt_pro8_trajout\WT_residue_pro9_rmsd.csv", mode="w", newline="") as file:
    writer = csv.writer(file)
    
    # Write the header row with residue IDs
    header = ["Frame"] + [f"Residue_{resid}" for resid in residue_rmsd.keys()]
    writer.writerow(header)
    
    # Write RMSD values for each frame
    num_frames = len(next(iter(residue_rmsd.values())))  # Get number of frames
    for frame in range(num_frames):
        row = [frame] + [residue_rmsd[resid][frame] for resid in residue_rmsd.keys()]
        writer.writerow(row)

print("RMSD values have been saved to 'WT_residue_rmsd.csv'")

----------------------------------------------------------------------------------------------------------------

import MDAnalysis as mda
from MDAnalysis.analysis import rms
import numpy as np

def extract_residue_rmsd(pdb_file, dcd_file, residue_id):
    # Load the structure and trajectory
    u = mda.Universe(pdb_file, dcd_file)
    
    # Select the residue of interest
    residue = u.select_atoms(f"resid {residue_id}")
    
    # Set the reference structure (first frame)
    reference = residue.positions
    
    # Initialize an empty list to store RMSD values
    rmsd_values = []
    
    # Loop through the trajectory
    for ts in u.trajectory:
        rmsd_value = rms.rmsd(residue.positions, reference)
        rmsd_values.append(rmsd_value)
    
    return np.array(rmsd_values)

----------------------------------------------------------------------------------------------------------------

import MDAnalysis as mda
from MDAnalysis.analysis import rms
import numpy as np

def extract_domain_rmsd(pdb_file, dcd_file, residue_ranges):
    """
    Extracts the RMSD of multiple residue ranges (a domain) over time.

    Parameters:
    - pdb_file: str, the path to the PDB file.
    - dcd_file: str, the path to the DCD file.
    - residue_ranges: str, the selection of ranges (e.g., "resid 100:200 or resid 250:375").

    Returns:
    - rmsd_values: np.array, the RMSD values for the domain over the trajectory.
    """
    # Load the structure and trajectory
    u = mda.Universe(pdb_file, dcd_file)
    
    # Select the group of residues (the domain)
    domain = u.select_atoms(f"{residue_ranges}")
    
    # Set the reference structure (first frame of the selection)
    reference = domain.positions.copy()
    
    # Initialize an empty list to store RMSD values
    rmsd_values = []
    
    # Loop through the trajectory and calculate RMSD for the domain
    for ts in u.trajectory:
        rmsd_value = rms.rmsd(domain.positions, reference)
        rmsd_values.append(rmsd_value)
    
    return np.array(rmsd_values)

----------------------------------------------------------------------------------------------------------------
# Each time you change the residue ensure you run this part
#  Specify your files for each simulation
simulations = [
    (r"r375k_pro9_trajout/traj.pdb", r"r375k_pro9_trajout/traj.dcd"),  #dark blue
    (r"r375l_v2_pro9_trajout/traj.pdb", r"r375l_v2_pro9_trajout/traj.dcd"),  #darkgreen
    (r"rw_lf_pro8_trajout/traj.pdb", r"rw_lf_pro8_trajout/traj.dcd"), 
    (r"w376a_v2_pro9_trajout/traj.pdb", r"w376a_v2_pro9_trajout/traj.dcd"),  #lightblue
    (r"wt_pro8_trajout/traj.pdb", r"wt_pro8_trajout/traj.dcd"),        #black
    # Add as many as needed
]

residue_id = "672"  # Specify the residue ID of interest, this can be a range of residues - specified by a colon, eg. 10:20 or several unlinked residues e.g. 100 110 114
#Can also specify two separate ranges at once: like "resid 100:200 or resid 250:375"
# Store RMSD values for each simulation
rmsd_data = {}

for i, (pdb, dcd) in enumerate(simulations):
    rmsd_values = extract_residue_rmsd(pdb, dcd, residue_id)
    rmsd_data[f"Simulation_{i+1}"] = rmsd_values
----------------------------------------------------------------------------------------------------------------
#or if you are trying to run for several residues in two separate ranges (e.g. 120:150 and 180:200)...

# Specify your files for each simulation
simulations = [
    (r"R375K/traj.pdb", r"R375K/traj.dcd"),   #dark blue
    (r"R375L/traj.pdb", r"R375L/traj.dcd"),   #darkgreen
    #(r"RW_LF/traj.pdb", r"RW_LF/traj.dcd"),
    (r"W376A/traj.pdb", r"W376A/traj.dcd"),   #lightblue
    (r"WT/traj.pdb", r"WT/traj.dcd"),         #black
    # Add as many as needed
]

residue_ranges = "resid 1:89 or resid 300:337"  # Specify the domain as two residue ranges

# Store RMSD values for each simulation
rmsd_data = {}

for i, (pdb, dcd) in enumerate(simulations):
    rmsd_values = extract_domain_rmsd(pdb, dcd, residue_ranges)
    rmsd_data[f"Simulation_{i+1}"] = rmsd_values

----------------------------------------------------------------------------------------------------------------


import matplotlib.pyplot as plt
import numpy as np
import os 


# Define the colors you want to use for each simulation
colors = ['darkblue', 'darkgreen', 'lightsteelblue', 'orange', 'black']      # Add more colors as needed

# Plotting the RMSD for each simulation with specified colors
plt.figure(figsize=(8, 5))

for i, (pdb_file, dcd_file) in enumerate(simulations):
    # Extract the directory name (or file name) to use in the legend
    sim_name = os.path.basename(os.path.dirname(dcd_file))  # Use directory name
    color = colors[i % len(colors)]  # Cycle through the colors if needed

    # Extract the RMSD values for the specified residue
    rmsd_values = extract_residue_rmsd(pdb_file, dcd_file, residue_id)

    # Plot RMSD values for this simulation
    plt.plot(rmsd_values, label=sim_name, color=color)

# Set x-axis limits to avoid extra white space (no white space on both ends)
#plt.xlim(0, len(rmsd_values) - 1) # Adjust based on the number of frames
plt.xlim(0, 510)
plt.ylim(0, 35)

# Add labels and title
plt.xlabel("Frame", fontsize=20, labelpad=20)
plt.ylabel("RMSD (Å)", fontsize=20, labelpad=20)
plt.title(f"Chain B F377", fontsize=24, pad=20)

# Increase font size for the ticks on both axes
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)

# Display the legend with the correct simulation names
plt.legend(loc="best", fontsize=12)

# Enable the grid for better visualization
#plt.grid(True)

# Save the plot as a PNG file with 500 DPI resolution
plt.savefig(' ChainB_F377_plot_4.png', dpi=500)

# Show the plot
plt.show()
